#!/usr/bin/runhugs
--created on 2009-03-18
--Playing with the normal statistical distribution.

--import Test.QuickCheck

sqr x = x*x

--approximation by k-th degree Taylor polynomial
erf_ z k = 2 / sqrt pi * summa where
  summa  = sum [z/(2*fromIntegral n+1) * pr n | n <- [0..k]]
  pr n   = product [-sqr z/fromIntegral i     | i <- [1..n]]

--non-continuous hack!
--monotonic increasing function, @-inf=-1, @0=0, @inf=1
erf z | z < -4 = -1
erf z | z >  4 =  1
erf z          = erf_ z 100

--probability density function for standard normal dist.
pdf 0 1 x = exp (-sqr x/2) / sqrt(2*pi)

--x<mu: mon. inc., x>mu: mon. decr., @-inf=0, @inf=0,
--max(x: pdf mu sigma x) = mu
--forall mu, s1, s2: s1>s2 <=> pdf mu s1 mu < pdf mu s2 mu
--forall mu, sigma, x: pdf mu sigma x >= 0
--probability density function for normal distribution
pdf mu sigma x | sigma>0 =
  pdf 0 1 ((x-mu) / sigma) / sigma

--monotonic increasing function, @-inf=0, @inf=1
--forall mu, sigma: cdf mu sigma mu = 0.5
--forall mu, s1, s2, x<mu: s1<s2 <=> cdf mu s1 x < cdf mu s2 x
--cumulative distribution function for normal dist.
cdf mu sigma x | sigma>0 =
  (1 + erf ( (x-mu) / (sigma*sqrt(2)) ) )/2

--regression tests for this program
failedTests        = filter off results where
  off (_, y, fx)   = abs(y - fx) >= 0.00001
  results          = map eval $ zip [0..] cases
  eval (n, (y, f)) = (n, y, f())
  cases = [
    (4,       (\() -> sqr 2)),
    (0,       (\() -> erf 0)),
    (0.5205,  (\() -> erf 0.5)),
    (0.8427,  (\() -> erf 1)),
    (1,       (\() -> erf 9)),
    (-0.5205, (\() -> erf (-0.5))),
    (-1,      (\() -> erf (-9))),
    (0.5,     (\() -> cdf 3 (sqrt 0.2) 3)),
    (0.98733, (\() -> cdf 2 (sqrt 0.2) 3)),
    (0.84134, (\() -> cdf 100 15 115)),
    (0.97725, (\() -> cdf 100 15 130)),
    (0.89206, (\() -> pdf 3 (sqrt 0.2) 3)),
    (0.07322, (\() -> pdf 2 (sqrt 0.2) 3)),
    (0.02659, (\() -> pdf 100 15 100)),
    (0.01613, (\() -> pdf 100 15 115))
    ]

main = print failedTests

--prop_cdf05 :: Double -> Double -> Property
--prop_cdf05 x y = y/=0 ==> cdf x abs(y) x == 0.5
--main = quickCheck prop_cdf05

--prop_pdfs1s2 :: Double -> Double -> Double -> Property
--prop_pdfs1s2 mu s1 s2 = s2>0 && s1>s2 ==> pdf mu s1 mu < pdf mu s2 mu
--main = quickCheck prop_pdfs1s2
